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ABSTRACT 



In this paper we explore the maximum precision attainable in the location of 
a point source imaged by a pixel array detector in the presence of a background, 
as a function of the detector properties. For this we use a well-known result from 
parametric estimation theory, the so-called Cramer- Rao lower bound. We develop 
the expressions in the 1-dimensional case of a linear array detector in which the 
only unknown parameter is the source position. If the object is oversampled by 
the detector, analytical expressions can be obtained for the Cramer-Rao limit 
that can be readily used to estimate the limiting precision of an imaging sys- 
tem, and which are very useful for experimental (detector) design, observational 
planning, or performance estimation of data analysis software: In particular, we 
demonstrate that for background-dominated sources, the maximum astromet- 
ric precision goes as -B/-F^, where B is the background in one pixel, and F is 
the total flux of the source, while when the background is negligible, this preci- 
sion goes as F^^. We also explore the dependency of the astrometric precision 
on: (1) the size of the source (as imaged by the detector), (2) the pixel detec- 
tor size, and (3) the effect of source de-centering. Putting these results into 
context, the theoretical Cramer- Rao lower bound is compared to both ground- 
as well as spaced-based astrometric results, indicating that current techniques 
approach this limit very closely. It is furthermore demonstrated that practical 
astrometric estimators like maximum likelihood or least-squares techniques can 
not formally reach the Cramer-Rao bound, but that they approach this limit in 
the 1-dimensional case very tightly, for a wide range of S/N of the source. Our 
results indicate that we have found in the Cramer-Rao lower variance bound a 
very powerful astrometric "benchmark" estimator concerning the maximum ex- 
pected positional precision for a point source, given a prescription for the source. 
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Introduction 



Astrometry relies on the precise determination of the relative location of, usually, point 
sources. The estimation of the precision with which these measurements can be done, both 
from an empirical, as well as from a theoretical point of view, has been the subject of 
various pape rs. Seminal work, as appli ed to stellar images recorded on ph otographic plates. 



are those of 



van Altena fc Auer (1975) and 



refinements by 



Lee and van Altena 



Auer and van Altenal (119781 ). with further 



(119831 ). in which statistical estimations for the precision 
of the position of stellar images were compared to the results from the actual fitting of 
stellar profiles measured using microdensitometer scans, through a classical least squares 
minimization technique assuming a Gaussian noise on the measured intensities. 



Nowadays, discrete digital detectors, such as Charged Coupled Devices (CCDs, iHowell 



(120061 )). being highly efficient area detectors , are widely us ed in astronomy for photometric 



astrometric and spe ctrosco p ic obs e rvations (IMackayl ul9St ), for t 



astrometry see e.g 



Monet 



(Il992[ ). 



LindegrenI ( 120051 ). and 



Howell 



l e spec ific use of CCDs in 
(120 131 )). This prompted 



KingI ( 119831 ) to carry out a similar analysis for CCDs, specifically for HST data, starting 
also from the assumption that a least squares minimization approach provides the best 
estimation of the relevant parameters. 



The studies by iLee and van Altenal (119831 ) and iKing] (119831 ) (see also IStond ( 119891 )) 



provide estimates of the statistical uncertainties of the fitted parameters, given a noise 
model for the data. However, a related question, less often addressed, is what would be 
the maximum attainable precision with which one could expect to estimate the astrometric 
position of a source, given a prescription of the detection process. This question constitutes 
a central aspect to astrometric work. For example, in situations, when the detector nearly 
critically samples the light distribution for point sources rendered by the telescope optics 
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S) 



(case, e.g., of the HST-WFPC imageio), the question may arise on how well in principle the 
flux and the position of a point source located on some unknown background may be jointly 
estimated. The answer may be needed in instrumental design, for planning observation 
campaigns, or for checking the quality of data analysis algorithms. 

In this paper, we concentrate, precisely, on deriving a lower bound to the expected 
astrometric error for the position of a source in one dimension (hereafter referred to as 
1-D), as specified by the variance of the position itself, for the kind of data expected 
from astronomical C CD-observa t ions. Sg r ne sem inal work in this area, u sing the so - called 



Cramer- Rao bound ( 



akobsen et al. 



iao J 



(11992h . 



19451 



Cramer ( 



Zaccheo et al. 



(120041 ). More recently. 



Lindegren 



mm 



19461) )■ h a s been presented 



Adorj f ll996h . 



by 



Lindegren 



Lindegren! ( 



fll997h and 



1978h . 



Bastia: 



iani 



( I2OIOI ) has published a review paper exploring the 



astrometric bounds, using some of the Cramer-Rao prescriptions, but focusing mostly on 
space-based near-diffraction limited imaging. A parti c ularly relevant and inspirational 



Winickl (119861 ) which was focused, however. 



early work on the subject is that contained in 
on read-out noise limited devices and applications. We use Winick's paper as a starting 
point of our discussion, expanding it to sky-background limited astronomical observations, 
and develop it further to explore the general Cramer-Rao bound for astrometry in 
the 1-D case. We also elaborate on some useful approximations that explicitly expose 
the dependency of the expected minimum astrometric uncertainty as a function of the 
parameters characterizing the source and the detector. Some interesting open astrophysical 
problems that require high - accura cy astrometric measurements have been recently described 



in the book by ivan Altenal (120131 ). One may use the maximum expected accuracy estimates 



developed in this paper to determine whether or not a particular astrophysical problem 



^Instrument handbook available at http://www.stsci.edu/hst/wfpc2/documents/handbook/IHB_17.htn: 
last accessed on December 2012. 
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- and its associated accuracy requirements - may or may not be tackled with certain 
instrumental set-up and observing conditions. 

2. Preliminaries 

To give a formal context, in this section we introduce the basic setting of parameter 
estimation, as well as some related concepts and definitions, that will be used throughout 
the paper. 

2.1. Parameter estimation and the Cramer-Rao minimum variance bound 

The parameter estimation problem at hand can be presented in one-dimension, in 
general terms, as follows: Let us consider a collection of Jj (with i = l...n) independent and 
identically distributed realizations of a random variable. In this setting, it will be assumed 
that the {!{ : i = l,...,n} measurements follow an underlying probability (density - if 
continuous, or mass - if countable) function, denoted by fg, which depends upon a certain 
target (unknown) parameter 6. Then, the parameter estimation problem reduces to find a 
prescription (or statistics) such that the function 6{Ii, ..., /„) is a good approximation of the 
underlying parameter 9 that generated the data. 

A standard criterion adopted in statistics to estimate 6 is to consider the rule of 
minimum variance (denoted by Var), given by: 



^0 = argmin\/ar(^(/i,. ..,/„)) 
e 

[•gminE7,,...,7„^/n fe{h, ...,/„) - 6 
v 



ars 
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where E is the expectation value of the argument and "argmin" represents the 
argument that minimizes the expression. Note that in the last equality we have assumed 
that ^0 is an unbiased estimator of the parameter (i.e., that K{6) = 6), so that under this 
rule we are implicitly minimizing the mean square error of the estimate with respect to the 
hidden true parameter 6. 

Unfortunately, the general solution of equation ([1]) is intractable, as in principle it 
requires the knowledge of 6, which is the essence of the inference problem. However, 
there are performance bounds that characterize how far can we be from the theoretical 
solution in equation ([T]), and even scenarios where the optimal solution can be achieved in 
a closed-form. One o f the most signif i cant r esults in this field is the Cramer- Rao minimum 



variance bound (JRad (119451 ) . ICrameii (Il946l )) explained below. 



Let 6{) be an unbiased estimator of 6. If we define the function L{Ii, ..., /„; 6) as the 
likelihood of the observation given the model parameter 6, and we have n independent 
random variables /« driven by the probability function fe, then the Cramer- Rao bound 
states that: 



varieih, ..., /„)) ^ Ej,_j^^fn e{h, ...,i„)-e] > 



If>(n) 



(2) 



provided that we satisfy the constraint: 



E/„...,/„./^. (^lnL(Ji,...,J„;^)l=0 



and where: 



(3) 



Io(n) = E 



h,-Jn-^fS 



dd 



lnL(/i,...,/„;i 



(4) 



is called the Fisher information of the data about 6. 

A powerful corollary of this result is that the minimum variance of any unbiased 
estimator that satisfies equation ([3]) is always going to be greater than the pre-specified 
quantity given by equation ([2]). 

Generally, this bound is not attained for the minimum variance estimator of 6, however 
there exists a necessary and sufficient condition that guarantees the existence of an estimate 
achieving the Cramer-Rao bound, Xg{n)^^. More precisely, if we can write (see, e.g.. 



Stuart, et al. 



bm pp. 12)): 

d\nL{Ii,...,Iri;0) 
dO 



A{d)-{e{h,...jn))-e) (5) 



then, it is certain that Kar(0(Ji, ..., /„)) = 1/X5)(n), as long as A{9) is a function 
exclusively of the parameter 6 (and, in particular it does not depend on the observables, 
/j)g. However, we must keep in mind that, unless the condition in equation ([5]) is satisfied, 
the minimum variance solution from equation ([T]) will have, in general, a variance strictly 
greater than Xg{n)^^. This important result will be further used in Section [51 

2.2. Position estimation: Astrometry 

The position estimation in astrometry is a slight variation of the classical parameter 
estimation problem introduced in Section 12. 1[ In this paper we focus on the 1-D version, as 
it is more easily tractable from the numerical and analytical point of view, while capturing 
all the key elements of the problem. The extension to the 2-D case will be dealt with in a 
forthcoming paper. However, as we shall see (see Sections 14.1.51 and [5!l. it is expected that 



^Furthermore, in this scenario, it can be easily shown that A{6) = Igln). 
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some generalizations are possible from our 1-D results which are likely to be approximately 
valid on the 2-D scenario in some cases. 

In the 1-D case we have an array detector with n pixels, in which we measure the fluxes 
{/j} per pixel. We will assume that the total expected (as opposed to measured) flux at 
each pixel on the detector, given by a function Aj(xc), will explicitly depend on the location 
of the source on the array, denoted by Xc-, which is the parameter we want to determine 
(equivalent to the unknown parameter 9 of Section [271]) . Of course, this flux is not measured 
directly, because the actual observations, Jj, are subject to noise. However, on photon 
counting devices, such as CCDs, the measured flux follow a Poisson noise distribution, i.e., 
the li are random variables driven by a Poisson distribution (this determines the probability 
mass function fo introduced in Section |2]), with expectation value given by Aj(xc). At 
this point, we note an imp ortant difference in approach to that adopted in the work by 
Lee and van Altenal (Il983[ ). in which they have assumed a Gaussian noise per pixel, valid 
for an analog detector, such as photographic plates. As we shall see, when the noise is 
Gaussian, a maximum likelihood parameter dete rmina t ion re duces to least squares, which 



Kind (119831 ) also adopted a least squares 



is what they indeed adopted. Note however that 

minimization, although for CCDs the noise is not Gaussian, and therefore a maximum 

likelihood solution is not equivalent to a least squares minimization (see Section [5]). 

In order to estimate the Cramer-Rao bound in this situation, we flrst need to verify 
that our likelihood function satisfles equation (j3]). The likelihood function of the 1-D array 
observations will be given by: 

L(/i, ..., /„; Xc) = fx,(x,)ih) ■ f\2(x,){h) ■ ■ ■ fx^{x,)iln) (6) 

where fx{I) = ^ p^ , since the Jj follow a Poisson mass function distribution with 
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i 



mean Xi{xc) q Then, we see that: 



~ ^dx^' "'^' = — f ^(/i-lnAi(x,)-Ai(a;c)-ln/i 



(7) 



dXi(Xr) sr^ dXAxr 



_ sr^ J 1 dXi{xc) Y^ 

. T ^i\Xc) (IX c . _, (IX c 

and, we indeed verify that E/^_ j^ f dx "' "''^ ) ~ ^ because E(Ji) = Ai(xc). Hence, 

we can apply equations (|2]) and (|4]) to obtain the following result: 



1 



Var{xXh. ..., In)) > :^TrT = 7— TT (9) 



E 



dxc 



^ XiiXc) 

For completeness, the derivation of the Fisher information about Xc, Xa,^(n), is presented 
in Appendix El 

In the following section, we provide a detailed analysis of this expression and its 
practical implications on astrometry. 

3. The astrometric Cramer-Rao minimum variance bound in 1-D 

In very general terms, the observed fluxes {/j} will have contributions from the source 
itself, as well as from a background. Correspondingly, the expected flux (in one pixel) from 
the source (which explicitly depends on Xc) will be characterized by a function Fi{xc), 



^Note that this estimation setting is different from the classical setting of Section [2Tl since 
we have random independent, although not identically distributed, samples. Nevertheless, 
it is simple to prove that equations (E]) to (j4]) still hold under this more general setting. 
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representing the flux (in photo-e~ on the detector) at pixel i, whereas the expected generic 
background will be denoted by Bi, representing the total (integrated) background (also in 
units of e^) at pixel i, and which includes contributions from the detector (read-out noise 
and dark-current, if any), and the sky background. We will assume that Bi does not depend 
on Xc- The total expected flux will thus be given by Xi{xc) = Fi{xc) + Bi. If we replace this 
expression for Aj(xc) into equation ([H]), we see that: 



Var{x,) > a^, = — ^^ (10) 

^— \ y dxc ^ ' 

i=i (Fi(x,) + Bi 



At this point, it is convenient to define a dimensionless, normalized, function gi{xc) 
such that Fi{xc) = F ■ gi{xc), where F is the total flux of the object (which is invariant to 
the actual value of Xc). In this case, the RHS of equation (ITOll can be written as: 



< = 7 rr- (11) 

" (^fe(^c)) 

i=i [F gi{xc) + Bi 



Winickl (Il98fih 



Note that this expression is similar to the Cramer-Rao bound derived by 
(his equation (35)). 

The function gi{xc) is determined by the Point Spread Function (PSF), which describes 
the distribution of (source) flux on the detector or, equivalently, the image profile across 
pixels (represented by the function $(a;)), integrated over pixel i (of width Ax) of the array, 
i.e.: 



gi{xc) = I <l>{x)dx (12) 

Ax 
2 
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Note that the PSF function $fx) is also normahzed, i.e. 



/+00 
<^{x)dx = l (13) 

-oo 



As long as the array length samples a significant fraction of the PSF, we can indeed 
identify F as the total flux of the star, since: 






J2Fi = F^gi{xc) = F^ ^{x)dx^F ^{x)dx = F (14) 

For practical purposes, since the detector array length greatly exceeds the PSF extent, 
this means that the source must not be too close to the array boundaries for this equation 
to be valid. 



3.1. Interpretation of the structure of the Cramer-Rao bound 

An equ ivalent to equation (llOl) in 2-D has been presented, in a slightly different 



manner, by 



Lee and van Altenal (11983) in the context of the expected astrometric 



accuracy on photographic plates (their equation (9)). Indeed, it is easy to see that their 
A?- = Fij{xc, Vc) + Bij = cr% + cr| , where we have assumed that both, the source and the 
background, follow Poisson noise (in the case of CCD detectors this is valid when the fluxes 
are measured in photo-e^), and where a\ and cr| are the variances in the source and 
the background, respectively, at pixel z,j. We emphasize however that their equation (9) 
has been derived in the case of a (weighted) least squares solution (their equation (5)). 
This equation is valid for an analog detector, such as photographic plates, where the 
p hotographic densities i n eac h pixel are assumed to follow Gaussian noise (see equation (4) 



m 



Lee and van Altenal (1l983l )). However, as it will be demonstrated in Section [5]), in the 
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case of digital detectors, where the noise follows a Poisson distribution, no parameter 
estimator can formally reach the Cramer-Rao bound. 

We note that, for Poisson noise, the standard deviation of the signal is the square root 
of the signal and, thus, the variance is the signal itself. Thus, the interpretation of the term 
^ii^c) + Bi in equation (TTUj) as the variance of the counts (in e~) is important, as indicated 
in what follows. It is often more convenient for evaluation purposes to express F and B in 
terms of "counts" on the detector (from now on referred to as "ADUs" - Analog to Digital 
Units), rather than in e~, by introducing the so-called (inverse-)gain of the detector G 



in units of e /ADU (see, e.g. 



Gillilandl ( 119921 )). In this case, the source and background 



fluxes, F and B, are given hj G ■ F and hj G ■ B respectively, where F and B are in 
units of ADU. On the other hand if ap, is the rms deviation at pixel i (in e~) and a p. is 
the equivalent quantity in units of ADUs, then it is true that a p. = G ■ a p. , and similarly 
for the rms deviation on the background. By replacing these unit conversions into either 
equation flTOj) or f fTTj) we see that we can express those equations with the source and the 
background measured in units of either e~ or ADUs, in the sense that: 



^CK = r-- z^ (15) 



n 

>: 

i=l 


(S(-c))' 


(4, + 4) 

1 


n 


(g(-=))' 



(16) 

u «. + o 

To evaluate the Cramer-Rao bound if we have empirical measurements of fluxes and 
variances on a detector (ADUs), it would obviously be more convenient to use equation (iT6l) . 
However, when performing numerical experiments for a given detector set-up (as done in 
this paper. Section l531) . where we only specify flux levels for the source and the background. 
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and where we calculate the variances associated to them, we would instead need to use 
equation f llSp (see also equation f l^T]) ). This is because we know that if the flux is in 
expressed e~, then the variances are equal to the flux (but this is not the case if the 
flux is measured in ADUs). We further note that equations ( !T5|) or ( IT6|) suggest that the 
Cramer- Rao bound represents a mean "uncertainty over the derivate of the signal" , since: 





^ca = — TT = TT^ (17) 



where the <> stands for a classical type of harmonic mean over the pixels, and where 
0"? = cr| + 0"! . This provides an interesting connection with what will be presented 
in Section [5.31 We note that if we have a function of one variable, y = f{x), thero 
ffy ~ (^) ■ a^. If we have n independent measurements, indicated by {xi,yi), each with 
uncertainty (cr^., ay-) we know from the propagation of errors in a least squares sense that the 
minimum variance for the weighted mean x would be given by cr| — ^ 



X " 1 



(see, e.g.. 



1=1 ^' i=i 




Meyer! ( 1l992l ). Chapter 10), which is equivalent to the Cramer-Rao bound if we 
identify fi = Aj and, since the errors follow a Poisson distribution, then cr^. = Aj. We 
emphasize that this analogy is valid only in the 1-D case (otherwise we would need to 
include the partial derivatives of /, and the variances in all the parameters). Indeed, as we 
shall see in Section El the least squares does not reach, in general, the Cramer-Rao lower 
bound. 

An interesting aspect of equations ( ITOj) or (ITTI) is that the positional uncertainty is 



"^Using a Taylor expansion around a point and assuming that / is sufficiently regular or 
smooth around that point. 
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going to be dominated by the region near the center of the PSF, where its derivative is 
steeper. Including regions far from the central core of the PSF, where ^(a^c) ~ 0, will not 
contribute to a decrease in the uncertainty and, instead, will only deteriorate the overall 
S/N of the source by incrementally adding more noise than signal (see the paragraph 
following equation (J28l) ). 



3.2. The Cramer-Rao bound for a Gaussian source 



It is evident that very little progress can be made in the estimation of the Cramer- Rao 



bound from equation ( ITTj) . unless we specify a shape for the PSF. Various anal ytica. 



have been proposed for the PSF of a p oint source as i maged by grouii d -base d (jKing 



forms 



and space-based detectors fJKing 



19831 ) (but, see also 



Bendinelli et al 



197ll) 



(|l988|)). Without 



loosing too much generality, in this study we will adopt a Gaussian function which seems to 
be a good representatio n of the PSF, at leas t from the stand-point of astrometric accuracy 



on ground-based data (JMendez et al. 



(J2010l )). and which allows some simple analytical 



manipulation (see Section WT\i . Under this assumption, we would have: 



^{x) 



'iTxa 



(x — Xc) 



arcsec 



:i8) 



where we adopt to measure x, Xc and a in units of arcsec. 

In this case, it is easy to show that the derivative of gi in equation flT2l) can be written 
in a closed-form, as follows: 



dgt 



[Xr. 



/27f(T 



-7(a;, ) 



-7(4) 



arcsec 



(19) 



where: 
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7(a;) 



I tXj tXj f> 



2a^ 



(20) 



with X- = Xi — ^ and xf = Xi + ^. 



Combining the results of equations ( fT2l) . ( ITSl) . ( fT9|) and ( !20l) in ( fTTj) and converting to 
ADUs, we finally arrive at the following exact expression for the Cramer-Rao lower-bound 
in 1-D for a Gaussian PSF: 



"^CR = 27ra^ ■ 



B 



GF^ 



E 



-7(a;, ) _ p-li^t) 



(21) 






Ax 
2 



-^(") dx 



Xi 



where we have assumed, for simplicity, that the bac k groun d is uniform (and equal 



Winickl (119861 ). this expression makes it 



to B) under the PSF of the object. As noted by 

explicit that the Cramer- Rao bound depends on F and B separately, and not just on the 

ratio F/B. If we adopt a and Ax in unit of arcsec in the sky, then the square-root of 



equation (1211) gives us the Cramer- Rao bound in units of arcsec directly. 



In Figure [T] we show th e resu 
setting proposed by 



ts of evaluating equation (12T1) under the experimental 
Winickl (Il986l ). i.e., assuming a fixed set of F and B values (and 
therefore a constant ratio F/B), for different values of the detector pixel size Ax. In this 
figure we introduce the "Full- Width at Half-Maximum" (FWHM), usually termed as 
"image quality" at astronomical observing sites, which is rel ated to t he Ga ussian a through 



FWHM = 2a/21ii2 a. Figure [His equivalent to Figure 1 in 



Winickl mm . 



We note that, for a given continuous pixel x-coordinate on the array, the corresponding 
(integer) pixel ID, i, on the array, is given by: 
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i = INT {x + 0.5) (22) 

where the function INT represents the integer part of the argument. In this paper we 
adopt that pixel ID z = 1 has pixel coordinates 0.5 < x < 1.5, pixel ID i = 2 has pixel 
coordinates 1.5 < a; < 2.5, and so on, following the convention of the IRAF packagqj. In this 
scheme, each pixel has width 1.0 in pixel x-coordinates, centered at x = FLOAT (i), with 
upper/lower pixel boundaries given by x± = FLOAT (i) ± 0.5 (where the FLOAT function 
converts an integer into a real number). The relationship between pixel x-coordinates and 
"physical" coordinates (projected onto the sky, as measured from the origin of the array), 
is given by px = [{x — 0.5) • Ax] arcsec. 

3.3. The Cramer-Rao bound and dithering in undersampled images 

An interesting feature is the effect on the predicted Cramer-Rao bound of pixel 
de-centering of the source. Figure [1] shows the effect of pixel de-centering for a 
FWHM = 0.5 arcsec, as a function of Ax. The fact that the Cramer-Rao b ound depends 
on the location of the source itself was already pointed-out by 



Winickl fll986f ). but this is 



evident from equations (TTOl) . (ITT]) or (pTj) . all of which explicitly depend on Xc- Of course, 
the effect is symmetrical with respect to the pixel center, so the impact of a de-centering 
of, e.g., +0.125 pix with respect to the pixel center is exactly the same as that of a 
de-centering of -0.125 pix, and (as long as the array properly samples the source), i.e. the 
effect is periodic (such that the Cramer-Rao bound is the same if the source is placed at 



^IRAF is distributed by the National Optical Astronomy Observatories, which are oper- 
ated by the Association of Universities for Research in Astronomy, Inc., under cooperative 



agreement with the National Science Foundation, see |http://iraf.noao.edu/ 
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X ± FLOAT [n), where n is an arbitrary integer) 



The important role of image de-centering in the cas e of undersampoled 



Anderson fc Kind (12000h . They 



images (such 



as those of HST, see below), has been demonstrated by 
refer to the "pixel— phase error" as the systematic error in the derived center if a centering 
algorithm is used that does not properly account for the distribution of signal among pixels 
when the source is not centered in a pixel. On the other hand, the Cramer-Rao bound 
described here represents instead the precision, or random error, statistically attainable 
with an ideal, unbiased, centering algorithm. 

The dotted line on Figure [1] is centered on a pixel, near the center of the array, the 
double-dotted-dashed line is for an offset of 0.125 pix, and the dot-dashed line is for a pixel 
de-centering of 0.25 pix. As it can be seen, the loss of astrometric accuracy as the pixel size 
increases, is less severe when the target is not at the center of a given pixel, but, rather, 
when it is offset from it. This is actually an intuitive result: For a given pixel size, when 
the source is not centered on a pixel, its flux is spread among more neighboring pixels, 
and therefore the source can be located more precisely. This result also implies that, for 
under-sampled systems, it is a good practice to "dither" the source a bit (even a fraction of 
a pixel) so that, in the end, the average Cramer- Rao bound is better than that if the source 



were located at the cente r of a pixel, a well 



flAnderson fc Kind fl2000h . 



ingwn technique applied, e.g., to HST images 



Fruchter fc Hook 



(I2OO2I )). To quantify the effect of dithering. 



in Table [T] we show the Cramer-Rao bound for a source with a FWHM = 0.5 arcsec 
observed through detectors of pixel size Ax = 0.5, 0.7 and 0.9 arcsec, and when the source 
is centered and slightly de-centered by the amounts indicated in the table. As it can be seen 
from this table, for a pixel size that matches the FWHM, the change on the Cramer-Rao 
bound is small as a function of pixel offset. In the case of a 0.7 arcsec pixel size, a dither 
pattern including offsets of 0.125 and 0.25 pix, plus a central pointing, yields an average 
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Cramer-Rao bound of ~6.0 mas, which represents an almost 18% improvement over the 
(single) centered Cramer-Rao bound. For a 0.9 arcsec detector the effect is even more 
dramatic, yielding an almost 40% improvement. 

Figure [1] also shows that, while the dithering technique offers a better asymptotic trend 
among all the values at large Ax (i.e., in the low-resolution regime), eventually, when Ax 
becomes too large, the astrometric accuracy deteriorates regardless of the relative position 
of the source with respect to the center of a pixel (even with dithering), as expected. 

4. Astronomical application of the Cramer-Rao bound 

We must note that F in equation (l?I]) is the total source flux, which is independent 
of the pixel size on the detector, whereas B is, instead, the background level in one pixel. 
Therefore, as Ax becomes smaller and smaller, the total contribution from the background 
under the PSF of the source increases steadily (because B is fixed, and the number of pixels 
under the PSF increases), and the positional precision deteriorates (for a specific connection 
with the S/N ratio of the source, see equation fl28|) and the comments that follow that 
equation). On the other extreme, as Ax increases, we loose resolution and the positional 
precision also deteriorates. We get a "valley" which determines an optimum region Ax for 
a given set of F, B (and a). While this setting is interesting in certain applications where 
the value of B is independent of the pixel size (e.g., military or day-time applications where 
the readout of the array is very fast, and the background is dominated by the electronics 
of the device, or when we are dominated by dark-current, see equation (!23|) ). the situation 
in astronomical applications is quite different: In this case, the long readout times imply 
that, in most cases, the background is dominated by diffuse light coming from the sky, and 



not from the detector, and in this case the bac kgrounc 



of the pixel size, as assumed in the analysis by 



i n a g iven pixel is not independent 



Winickl ( 119861 ). The setting for evaluating 
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equation fl2Tl) must then be adapted to our case of interest, this is done in the next section. 

As indicated in the previous paragraph, the correct expression for the (constant as a 
function of x, see equation (pTj) ) background B contained in one pixel, is given in this case 
by: 



fl.;,A.^£±|2^[ADU] 



(23) 



where fs is the sky background (in units of AD Us/arcsec), w hile D and RON are 



the dark-current and read-o ut noise of the detector (IHowell 
units of e~. In the paper by 



20131 pp. 222)a, per pixel, in 



Winickl ( 119861 ) it was assumed that /s ~ (true for very short 
exposure times), and we can see that, indeed, in this case, the background is independent 
of Ax. In what follow we will neglect the contribution from dark-current, which in current 
CCD detectors is negligible. 



With this new prescription for B, and in order to evaluate the RHS of expression (!2Ti 
for some astronomically interesting situations, it is is also worth to develop some easily 
measurable form of "signal" and "noise" for our Gaussian source, as observed through our 
CCD detector. In this case one could define the signal S, as: 



S = G-F- 



^{x) dx [e 



(24) 



Xl 



where xi and x„ are suitably chosen (but arbitrary) apertures that include an 
appreciable fraction of the total flux of the star (we can not actually measure from — oo to 
+00 with a real detector, nor we want to do that since, in this formulation, the background 



^See also, e.g., http://www.ucolick.org/~bolte/AY257/s_n.pdf, last accessed on Decem- 
ber 2012. 
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will add- up to infinity over that aperture as well, see equations ( 126|) and fl27j) ). For the 
case of the Gaussian function adopted here, equation flTSl) . it makes sense to perform an 
integration of the PSF that is symmetrical with respect to the center of the source, centered 
at Xc, in which case the signal can be written as: 



S = G-F-P{u+) [e-] (25) 

where P{u+) is the probability integral] evaluated at u+ = {xu — Xc) /V^cr, and where 

U •^^C '^C •*^l' 

The total noise, A^, has contributions from the read-out-noise of the detector, the noise 



from the sky, and the noise from the s ource itse l: 
statistics (in e~), such that (see, e.g. 



Gillilandl f ll992h ): 



, all o f which are assumed to follow Poisson 



N= JS + A/pix (Gfs Ax + RON^) [e" 



(26) 



where N„^^ is the number of pixels under the same region in which the signal S was 
sampled, i.e., in the interval [xi,x„], which is given by: 



K 



Xu-xi 2^/2au+ u+ FWHM 



pix 



Ax Ax 7^2 

Combining equations ( j25ll . (l26ll and ( 1271) we see that: 



Ax 



(27) 



S 

N 



P(u^) ■ F 



,Ua 



Pju+yP _|_ 2V2u+ a 



G 



G Ax 



(/.Ax + ^) 



''The probability integral is defined as P{u) = -^ JJ'e ^ dv 



-22- 

^ +^ 28 

P{u+yF I u+ FWHM f'f A™ I RON'^ \ 
G ^ Vbr2G Ax \Js^-^ ^ G ) 

For example, for P = 0.9 (aperture containing 90% of the total flux), we have u^ ~ 1.164 
(the true "physical" aperture on the detector would be ig x FWHM ^ 1.40 x FWHM), 
whereas for P = 0.99 then m+ ~ 1.822 (or ^ x FWHM ^ 2.19 x FWHM). 
Because u+ increases faster than P{u^) (which is bound to a maximum value of 1.0), 
we see that, for a given source, background, and detector, the S/N computed from 
equation ( l28l) decreases as u+ increases beyond the main core of the PSF. For example, for 
Ax = 0.2 arcsec, G = 2 e'/ADU, RON = 5 e", /, = 2 000 ADU/arcsec, F = 5000 ADU, 
and FWHM = 1.0 arcsec, then for P = 0.9, S/N ~ 74, whereas for P = 0.999 (for which 
M+ ~ 2.33), S/N ~ 68. As was explained before, we note that equation ( |2T]) does not 
depend directly on the S/N. 

Equation ( 128|) is interesting since it explicitly shows that, as Ax becomes smaller and 
smaller, the RON term starts to dominate over the sky background in its contribution to 
the total noise, the impact of which, on the Cramer- Rao bound, has already been mentioned 
in Subsection 13.21 However, when Ax increases, the sky background becomes the dominant 
source of background noise, and the total noise becomes independent of the array pixel 
size. Also, this equation clearly shows the classical result that, as an image becomes more 
spread (larger FWHM, or worse image quality) the S/N deteriorates, for a fixed total 
flux F, because of the larger contribution from the sky and the (larger number of) pixels 
underneath the aperture: As we shall see, the FWHM has a very relevant impact on the 
Cramer- Rao bound (see, e.g., equation ( 145|) ). 

Figure [2] shows the result of evaluating equation (12T1) under the assumption of a 
background given by equation fl23|) for a set of representative values. An interesting point 
here is that, at very small values of Ax we still see the "upturn" in the Cramer- Rao 
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lower bound seen in Figure [H but it has a much smaller effect. Of course, the reason 
for this upturn is the prevalence of the RON over the sky background indicated in the 
previous paragraph, when Ax becomes extremely small. As we shall see (equation fH^ ). 
the Cramer- Rao bound goes as Ax~^ for small S/N and small pixels, a feature clearly seen 
in Figure |2l Otherwise we see a broad region that exhibits a rather smooth and steady 
decrease in positional precision when Ax becomes larger and larger, and a rather steep 
increase when Ax increases beyond the FWHM. The overall effects of pixel de-centering 
are qualitatively similar to those already presented in Figure [H and are thus not repeated 
in this figure. For very large S/N , equation fH^ predicts that the Cramer-Rao bound 
becomes rather insensitive to Ax, which also coincides with the behavior in Figure |2l 

An interesting prediction of equation (12T]) is that high-resolution imaging in low- 
background, even for under-sampled images (e.g., HST), is better than imaging with 
larger aperture ground-based telescopes, not under-sampled, due to the worse FWHM 
and higher-background of the latter, a well-known fact by people doing astrometry with 
HST (provided, of course, that systematic effects are well understood, e.g., a particularly 
challenging situation with HST data is t he account of t ime-d ependent charge-transfer 



eff iciency correct i ons, 
or 



or details see, e.g.. 



Bristow et al. 



teOOSi ). especially their Figure 4, 



Bristow et al.l (J20061 ). especially their Figure 10). For example, for the same detector 
parameters as those adopted in Figure [21 and F = 10 000 (which for a Gaussian PSF leads 
to maximum flux in the central pixel of ~ 1 700 ADU (see Section W7\\ and equation dS])), 
fs = 3 000 ADU/pix, and FWHM = 0.45 arcsec, the Cramer-Rao bound is ~1.7 mas (with 
Ax = 0.08 arcsec). Thes e (source fc image) y alues are similar to thos e of the QSOs used i n 



the astrometric study by 



Mendez et al. 



(12010[ l (see their Table 1) and 



Mendez et al. 



(1201 ih 



w hich demonstra t ed a single- measurement astrometric precision of 1.5 mas (see Section 3.2 



m 



Mendez et al 



( l2010f )) with the NTT (3.5m aperture) telescope and SUSI2 imager. On 



the other hand, for HST with fs = 30 ADU/pix and FWHM = 0.15 arcsec, then the 
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Cramer-Rao bound is ~0.2 mas (in this case Ax = 0.1 arcsec), whereas iPiatek et al.l (|2002| ) 
reported a single-measurement precision of 0.25 mas (in our calculation of the Cramer- Rao 
bound for HST we have approximately taken into account the aperture difference between 
the NTT and HST, and the different exposure times for the same QSOs adopted in these 



two studies, from Table 1 in 



Piatek et al. 



f l2002h ). 



4.1. The Cramer- Rao bound in the small pixel (high resolution) 

approximation 



Under certain circumstances, the summation in the denominator of the RHS of 
equation ( ITOll can be approximated into an integral, which allows us to explore the behavior 
of the Cramer-Rao bound in a more explicit manner. Indeed, we see from equations (IT2l) 
and flT8|) that the application of the mean-value theorem when Ax/a ^ 1 implies that 
Fi = F ■ gi{xc) ~ F ■ $(xj) ■ Ax. In this case, for a Gaussian PSF, it is easy to show that: 



^(x,)^F^(x,) 



dx, 



dXr 



I Xi Xr 



a" 



■F 



(29) 



Replacing the RHS of equation fl2I?]) into the RHS of equation flTUl) we have: 



^2 4 



X^i=l \-^i -^c) 



F2 



(30) 



F,+B, 



Let us have a closer look at the (dimensionless) Fi and F^ terms in the RHS of 
equation (130!) . For our Gaussian function we will have: 



X^ + ^ 



^(x) dx 



27ra 



2 (x-xc^ 

e 2<T^ dx 



(31) 
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which, in the small pixel size approximation becomes: 



Fi ^ ^^ ■ e ^^^ ■ Ax (32) 

V27rcr 

= -f^iax ■ e 2a^ (33) 

where -Fmax is the (dimensionless) maximum flux, which will occur at a certain pixel j 
that satisfies the condition Xj — Ax/2 < Xc < Xj + Ax/2. 

Equation (!32l) prompts us to define a new function, which describes the distribution 
across pixels, in the small-pixel approximation, of the square of the flux that appears in 
equation (l30l) . as: 

Ff = F^- e ^^~ ■ Ax (34) 

where F"^ is a proper normalization factor (in units of arcsec"^, see below). 

Combining equations (15^ and ( IMIl in equation ( 15U]) . and considering that Ax is very- 
small, we have: 



a„„ = -^r- ■ lim 5 = ^^ ■ 5 (35) 

tr (F,+5,) ■/— F(x)+5(x)^ 



Taking into account the definition of -Fmax on equation (1331) . from equation ( IMI) we 
will have that imax — ^j ^ -^^ ' ^^) where we have used the fact that, since Ax is very- 
small, then Xj ~ Xc- Replacing this approximation and equation fl33|) into equation f p5|) we 
get for the Cramer-Rao bound in the small-pixel approximation: 
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2 


a^ 






Ax 




-^max 


Ax 


P2 

-^max 


/•+00 


{x- 

( ~ 

-^max 


, r, ix-Tc) 


- dx 


J 




-Xc) e 






■e 2-^ +B[x) 





(36) 



where we note that (the trivial definition of) integral / has units of arcsec^. 

This formula is only valid under the assumption of small pixels, so, in a general 
application, equation (ITO|) should be used instead, or equation ( 12T|) for a Gaussian PSF. 
However, the truly interesting aspect of this equation is that it can be explicitly evaluated 
in two extreme cases: When the detection is dominated by the source, and when it is 
dominated by the background. This is done in the next subsections. 



4.1.1. Weak source 
When the background dominates, then Fmax ^ B{x). In this case we would have: 



2 _ {^-xcf ,^ 2 - <^-^g)^ 



^max ■ e ^^^ + B{x 



"-"^' ^ " dx^ ^"-"-! ' " dx (37) 



If we assume an approximately constant background under the PSF of the target then: 



1 f^°° . .9 (^-^c)^ . .AF o^ 



1=^1 (x- Xc? e~^~^^ dx = ^^^-"-^ (38) 

Therefore, in this case, replacing / into equation 0361) and re-arranging terms, the 
Cramer-Rao bound becomes: 



2 2 5^ 

(T^„ = —^ ■ ^^ ■ Ax ■ a 



A F]' 



max 
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B 



V27rln2 G -Fmax 



Ax ■ FWHM 



(39) 



The (pedagogical) use of this equation is that it allows us to draw some basic 
conclusions regarding the expected positional accuracy in this regime. First, because 
-^max ^ B{x) we predict, in general, a rather large positional uncertainty, as expected. 



due to the low S/N of the source. Furtherm ore, the accuracy wi 



to -Fniax' i^ agreement with equation (10) of 
first line of equation fl39|) with equation (7) in 



Lee and van Altena 



1 impr ove proportionally 



Auer and van Altena] (119781 )). Furthermore 



f|l983[) (c ompare also the 



as intuitively expected, the accuracy deteriorates for a larger background, coarser pixel size, 
or for lower-quality images (or sites), but relatively slowly: Only as the square root of these 
parameters. 

4- 1-2. Strong source 

In this case, the signal from the source dominates over the background, i.e., 
-^max ^ B{x), and the approximation for I becomes: 



+00 



n2 _i^-^c) 

[X — Xr) e ^ 



Fmax ■ e ^^^ + B{x] 



dx 



max "'-oo 



+ C)0 



{x — Xc) e 2a^ dx (40) 



Evaluating the definite integral we have: 



27r- -^ 



a 



-fmax 
Replacing this value for / into equation (l36l) we end up with: 



(41) 
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CR 



2 11a 

(Tr,„ = ■ — ■ Ax ■ a 

2tx Fmax 

^ ^ Ax-FWHM (42) 



4-\/vr In 2 G -Fmax 

— 1/2 

We see that, in this regim e, the ultimate positional a ccuracy is proportional to -^niax' 



similarly to what was found by 



Lee and van Altenal (119831 ) (see their equation (13)). Again, 



the accuracy deteriorates slowly for coarser pixel size and for lesser-quality images, but in 
this case the background level, formally, plays no role in the expected accuracy. 

4- 1.3. Limiting cases as a function of total flux 

Equations ( 139|) and ( H2|) are a bit misleading because, by definition, Fmax depends 
itself on the adopted value for Ax and, therefore, needs to be evaluated in each particular 
case. However, in the small pixel approximation, one can find an approximate relationship 
between the total flux F (which is independent of Ax) and -Fmax- Indeed, assuming, as 
done before, that, Xj ~ Xc, then: 



P /'^J + ^ ( - )2 

i^max = F, = -^^- / e"^?^ dx (43) 

F Ax 



V2n o 
Replacing equation (jH]) into equations fl39|) and fH2|) we have: 



(44) 




0F B FWHM^ if F <^ R 

2(2 1n2)3/2 ■ GF2 ' Ax i-i- r <^ D 
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Interestingly, we now see that the positional accuracy should deteriorate linearly with 
the FWHM of the image for strong sources, but it will not otherwise depend on the array 
pixel size. However, for weak sources, the dependence on the FWHM is not only steeper 
but, also, the finer the pixel size, the larger the expected positional uncertainty. This 
latter result should also be intuitive since, when we are dominated by the background, 
the more pixels we have under the PSF, the larger the contribution of the background 
will be, to the point of significantly perturbing the final positional accuracy. Equivalently, 
note that the term B/Ax is the background per unit area (see equation fj23|) ): The larger 
this value becomes, we indeed expect a larger positional uncertainty, as shown by the first 
line on equa tion ( H5|) . The behavior depicted by equation fH5|) is similar to that found by 



King] (119831 ): For sources where the background dominates, his equation (23) predicts that 



the positional uncertainty goes as yB/F, whereas when the background is negli gible, his 



equati on (24) sh ows that the positional un c ertain ty increases as 1/vF. Similarly. iLindegren 



fll978h (see also 



van Altena and Fomalont 



(12013 



pp. 127)) finds that for a seeing-limited 



image with no background, the limiting astrometric precision goes as FWHM/{S/N), which 
is equivalent to equation ( H5]) for the case F >> B, when S/N ~ -\/F, see equation (l28l)). 
All these predictions coincide with those from our equation ( |45|) . Finally, and as already 
noted in Section HI the general trends seen in Figure |2] are well explained by equation ( l45l) . 

4.1.4- Range of use of the high resolution Cramer-Rao bound 

Besides their qualitative usefulness, what is the range of applicability of equations (139|) 
and ( 142|) . or ( 145|) ? For illustration purposes, in Table [2] we compare the values predicted by 
these equations, with the "exact" prediction from equation (I2ip . for some representative 
values of the parameters. From this Table we can see that, even for relatively large 
pixels, in comparison with the FWHM, these equations predict very reasonable values in 
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comparison with the "exact" value, as long as we respect the conditions for F/ B under 
which equations fl5^ and f l42p were derived. The second and third lines show that, indeed, 
as predicted by equation f HS]) . the Cramer- Rao bound does not depend on B. The third and 
fourth lines show that in the high-S'/iV regime the Cramer-Rao bound goes linearly with 
the FWHM, while a comparison of the first and fifth line demonstrate that in the low-S/N 
case the Cramer- Rao bound varies as the ratio FWHM^/Ax, as shown by equation (I45p . 
Finally, in the last two lines of the table, we show a case when equation (1451) fails miserably, 
i.e., for an intermediate S/N value. 



4.1.5. A simplified extension to the 2-D Case 



It can be intuitively argued from either equations fl39|) and f l42|) . or equation fH5|) . 
that in the 2-D case the numerical factors in front of these equations will be somewhat 



altered, but their basic dependence on F and B should be ba sically maintained, as a 



indicated by a comparison of o ur resu l ts to the 2-D results by 



Auer and van Altenal (119781 ) 



ready 



Lee and van Altenal ( 119831 ) and iKind (119831 ) noted in Sections liXTl WT72\ and l4X3l For 
example, if the source is relatively symmetrical, such that we can replace an x-y integration 
(similar to that of equation (15^ but in the two array dimensions x, y) by a polar-radial 
integration, we will end up with a Ar instead of a Ax, or even something like a/ Ax ■ Ay for 
a non-square pixel array. Therefore, even though our present results are based purely on a 
1-D array, they have a very interesting predictive power, a point to which we will return in 
Section 15.31 (see also Table H]). 
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5. Comparing the Cramer-Rao bound with the performance of practical 

estimators 

Let us remember that the necessary and sufficient condition for the existence of an 
estimator that achieves the Cramer-Rao bound is that the likelihood function can be 
decomposed as in equation (|5]). Unfortunately, equation ([8]) does not offer, in general, the 
normal form of equation (|5]) , and consequently no estimator achieves the Cramer- Rao lower 
bound. 

To further explore this, it is illustrative to consider the high-resolution regime of 
Section WA] in which case equation fl29|) holds. Replacing this into equation (|8]), and 
re-arranging some terms, we have: 



rflnL(Ji,..., J„;x,) _ ^^ / _^^ _ 1 J . (^^ _^^) (46) 



dxc ^ a^ \Fi + B. 



i=l 



Note that although the expression in equation ( H6ll is reminiscent of equation ([5]), the 
factor attributed to A{6) in equation (146|) . is a function of the data {Ii,Xi), and therefore 
does not fulfill the decomposition in equation (JS]). Consequently, even under the high 
resolution approximation, there is no estimator that achieves the Cramer-Rao bound. This 
situation supports and justifies the adoption of alternative criteria for position estimation, 
maximum likelihood and the classical least squares being two of the more commonly used 
approaches adopted. These are reviewed in the following subsections. 

5.1. Maximum Likelihood 

Given the fiux values Ji, ..., /„, the maximum-likelihood (ML) estimate of the position 
Xc is obtained through the following rule: 
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Xc^,^(Ji,...,4) = argmaxlnL(Ji,...,4;a;c) (47) 

where "argmax" represents the argument that maximizes the expression. Imposing 
the first order condition on this optimization problem, it reduces to satisfying the 
condition dx ' "'^ ~ ^^ ^^'^' consequently, we can work with the general expression 

in equation ((71). We note, from that equation, that the term in the brackets given by 
Yll=i ^i{^c) = -F + -B is independent of Xc, and consequently its derivative is zero. Therefore 
the ML condition becomes: 



dlnL(/i,..,/„;3:e) ^ A 1 d\,{xc) _ ]_\^ j H^c) / _ ^ n _ n (AfA 

where we have used the high- resolution approximation, equation fl29|) . In general, 
this condition does not offer a closed-form solution (note that the dependency of Fj on 
Xc could be quite complex), consequently numerical gradient- descent methods need to 
be implemented. Nevertheless, in the signal-dominated regime, when F.i ^ Sj, it is 
straightforward to show that equation (l48l) reduces to the classical moment solution, i.e.: 



n 

'I ■ Xi 



*o„, = ^, (49) 



1=1 



In general, the formal solution to equation (148|) in the high-resolution regime offers a 
simple relationship to implement a recursive algorithm that satisfies the ML estimate, given 
by: 



-33- 



XcML = ^ (50) 

'^Wi{xc) ■ h 



where the weights Wj(xc) = 



H^c) 



Fi(xc)+Bi- 

To conclude this subsection, we must mention that it is well-known that in the setting 
of independent and ide ntically distributed measurements, the ML estimate is asymptotically 



unbiased and efficient (IStuart. et al. 



2004 . Chap. 18). In other words, as the number of 
samples increases, the variance tends to the Cramer-Rao bound and consequently the ML is 
asymptotically optimal (in the sense of achieving the minimum variance bound). However, 
the astrometry setting is different as the measurements (fluxes) follow a Poisson distribution 
with parameters that are position dependent (see Section 12^21) . and consequently the ML is 
not necessarily efficient in this statistical sense. 



5.2. Least Squares 

Given the flux values Ji, ..., /„, the (weighted) least square (LS) estimate of the position 
Xc is given by the following decision-rule: 



/T T \ ■ Sr^ l-'i Ctj(^c)) /■r:-,\ 

Xc^g{h,...,In)=argmm2_^ j-i—. (51) 



i=l 



bi{xc) 



I 



being ai{xc) = IE(/j) = Aj(xc) and bi{xc) = Var{Ii) = Aj(xc) for all ijj The unweighted 
LS estimate assumes instead that 6j = 1, where no variance normalization is considered at 



^Note that the ML reduces to the LS estimate when the data follows a Gaussian distri- 
bution. 
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all. 

If we define E^sih, ■■■,In] Xc) = J2^=i i^i ~ -^jI^c)) , then the first-order condition over 
the unweighted LS estimate imphes that: 

Hence, using again the high-resolution regime (equation ( 129|) ). this condition imphes 



that: 



2_^—-{Ii- K) ■ [Xi - xj = (53) 



dxc ^^-^ (t"^ 



and, consequently, the problem reduces to solve the equation: 



y^.Fj ■ {li - \i{xc)) ■ Xi 

^c^s = ^ (54) 

^F,-(/i-A,(x,)) 



i=l 



Note that this expression is the counterpart of equation fl50|) for the ML estimate, 
emphasizing that these two techniques offer different estimates of the position in the 
parametric setting of astrometry. Again the solution of (l54l) involves the use of numerical 
methods, as, in general, no closed-form solution can be derived from it. 

5.3. Numerical results 

In this section we present some results from numerical experiments for the standard 
deviation of the astrometric position of a 1-D Gaussian source sampled by a linear detector. 
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The goal of these experiments is to compare the performance of the ML and LS solutions 
described in the previous subsections, to the theoretical Cramer-Rao bound, under different 
assumptions regarding the detector and the source. 

Basically, we start by adopting a set of values for the gain and read-out noise of the 
detector, as well as its pixel size and number of pixels. For a source with a Gaussian PSF, 
we specify its width, the maximum flux at the central pixel and its center location (xc). For 
the background, we adopt a certain (fixed) value per pixel. The total flux of the source (in 
ADUs) is obtained through equation ( H3l) by direct integration of the Gaussian PSF given 
the adopted values for Ax, x^ and a. We generate many possible "observations" for the 
same combination of parameters, using a random-number generator dri ven by a Po i sson 



Press et al. 



fll992h . 



distribution using the poidev, gammaln, and rani routines explained in 

We note that we transform all the ADUs (source and background) to units of e^ using the 

adopted gain, before randomizing the data. When g enerating t he da ta, we also consider the 



read-out noise and the "digitization noise" (see, e.g. 



Gillilandl ( 119921 )). On output we have 



an array of pixel positions (1 to n) and their corresponding fluxes. 

After generating a (large) number of simulations for a given set of parameters, we 
compute the value of Xc using a ML and a LS (weighted and unweighted) procedure. 
Because we are operating in the 1-D case, where we are estimating the object's position 
exclusively, all the other parameters were assumed to be known, and fed to the routine that 
searched for the value of Xc, using either equation P7|) in the case of ML or equation fl5T]) 
in the case of LS. To estimate the optimum number of simulations required to obtain a 
stable solution, we computed a very large number of simulations for our first set, and then 
calculated the mean of the recovered Xc and its standard deviation, ax^, as a function of the 
number of simulations. The optimum value for the number of simulations should render a 
set of values {xc, cr^J that do not depend appreciably (within their statistical uncertainty) 
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on adding more simulations. For our simulations, this number turned out to be ~250 
simulations per set. 

In Table [3] we present the standard deviation from the simulations, Ox^-, for a number of 
representative cases, as well as the calculated Cramer-Rao minimum variance bound based 
on equation ([21]), denoted by ctcr. In all cases we adopted an array size of 100 pixels that 
properly covers the PSF, a pixel size of 0.2 arcsec, a detector RON = 5 e~, a (fixed) sky 
background of 300 ADU per pixel, and a Gaussian source with a FWHM = 1 arcsec. 

To our surprise, the results on Table [3] indicate that the performance of the ML and 
LS estimators are almost equal to the Cramer-Rao lower bound. This is a remarkable 
result, and it validates the predictive power of the Cramer- Rao bound and its use as 
a benchmark indicator in astrometry. From the theoretical side, we formally showed 
(Section m equation f l46|) ) that the Cramer- Rao bound for astrometry can not be achieved. 
Our numerical results suggest then that, both, the ML and the LS methods (which are 
widely used in astrometry), are very efficient in the sense of asymptotically (with the 
number of measurements) approaching the Cramer-Rao bound. Proving this conjecture is 
an interesting topic of further work which will help us to explain the results obtained, and 
further consolidate the use the Cramer-Rao bound for performance analysis. A possible 
explanation, valid in the 1-D case, has been proposed in Section 13.11 This explanation 
however is not valid in a general situation, where one needs to simultaneously estimate 
astrometric and photometric parameters, the subject of which will be further explored in a 
forthcoming paper. We also note from Table E] that the ML is as good as, or even better in 
some cases, than the LS and WLS estimators. 

We finally compare the predictions of the Cramer- Rao lower bound to the performance 
of some 2-D digital centering algorithms applied to simulated stellar images, reported by 



Stond ( 119891 ). As argued in Section 14.1.51 we should expect that our results from the 1-D 
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case sho uld no t signi ficantly differ from those of a 2-D case when the source is symmetrical. 



Indeed, IStond (Il989[ ) adopted for his simulations symmetrical Gaussian sources (his 
equation (1)) on a fiat background. We have read approximately from his Figures 2 to 5 
the results for the minimum rms (over different profile fitting methods) of the positional 
uncertainty from his 2-D numerical simulations for some representative values of the total 
fiux F, this are tabulated in our Table S] under the label o"st, in units or arcsec. For each of 
these values we have computed the predicted Cramer- Rao bound u sing our equa tion (12T|). 



and adopting the same values for F, B, FWHM, and Ax used by IStond ( Il989l ) for each 
of his simulations. These values are denoted by o"cr in Table HI also in units or arcsec. 
As can be seen from this Table, our predictions are extremely encouraging: In all cases 
our computed Cramer-Rao bound is smaller than (although typically close to) the results 
from the "measured" standard deviation of the position (derived from the 2-D simulations). 
Furthermore, as indicated in Section |5TT] (specially equation f l50|) ). when F » B we would 
expect that the ML estimate approaches the moment solution, which is i ndeed the cas e as 



seen from Table H] (see also Figures 2 and 3 for large "counts in image" in IStond ( 1l989h ). 



Conclusions 



Our results indicate that we have found in the Cramer-Rao lower variance bound 
a very powerful astrometric "benchmark" estimator concerning the maximum expected 
positional precision for a point source, given a prescription for the source, the background, 
the detector characteristics, and the detection process. 

We regard as particularly interesting, pedagogical, and as a "back of an envelope" 
estimation tool, the set of equations (HSj) which, albeit derived in the high-resolution regime, 
are actually quite resilient to this condition, and thus provide reasonable expectations for 
the astrometric precision. 



-38- 

In a forthcoming paper we will formally extend our analysis to the 2-D case for the 
estimation of (x, y) coordinates in a more realistic pixel array. It has been argued however 
(see Sections 14.1.51 and I5.3I) that, as long as the PSF is reasonably symmetric, the results on 
equation (H5!) are not likely to change much in this case. In particular, the result that for 
background-dominated sources the Cramer-Rao lower bound goes as B/F'^, while when the 
background is negligible, this maximum achievable precision goes as F~^ should still hold 
in the 2-D case. 

Also, it would be interesting to study the sen sitivity of th e Cramer-Rao bound to other 
PSF shapes. We note however that the results by iKingl (jl983l ). which coincide with our own 



result s (see Section 14.1. Sp . have been derived from a different PSF (his equation (10)). 



King 



(119831 ) argues that, by using a Gaussian instead, only minor differences in the numerical 



coefficients in his equations (23) and (24) appear. 

Other factors to consider in a more realistic application include the existence of a 
strong gradient in the background under the object's PSF, which will tend to bias the 
derived astrometry, and possibly also affect the Cramer-Rao estimate (throughout this 
paper we have assumed that the background is uniform, but an extension to a highly 
variable background is straightforward to implement, starting from equation (TTTl) ). Also, in 
the case of severely under-sampled images, the issue of int ra-pixe l resp onse function has 



Adorj (ll996|)), and should be 



been demonstrated to be significant (see, e.g.. Figure 1 in 
included in the analysis. 

Eventually, one would like to be able to compute the Cramer-Rao bound in cases 
where there is a simultaneous joint estimation of the photometric (source and background), 
and astrometric (position, width of the PSF) parameters, which would for sure require 
non-linear, possibly iterative, numerical methods, in which careful attention should be given 
to numerical stability issues, and proper handling of observational errors - this constitutes 
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the long-term goal of our research. 
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A. Derivation of the Fisher information about Xc 

We start from the expression in equation ([7]): 



d\iiL{Ii,..Jn]Xc) _sr^ J 1 d\i{xc) ,.. 

dxc -^ * Ai(xc) dxc 

where, without any loss of generality, we have assumed that Yl^=i -^il^c) is constant, 
independent of the parameter x^ (see details on Section [3?T]) . Then, the Fisher information 
about Xc can be obtained as follows (see equation (jl])): 



1,. in) 



= E,,,...,,„(^(^lnL(/i,...,J„;x,))'^ 

m I Y^" ^n Ij-Ij d\i{xc) dXj{xc) 

— '^Il,...,In \yl^i = l l^j = l X^(^Xc)-\j{Xc) ' dXc ' dxc 

W I V f ^i d\i{xc) \ ) I m (sr^ V ^i'-^j dXjjXc) d\j{xc) \ 

^Il,■■■,h^ \l^i \^\^[X,) ' dXc J j ^ ^h,-,Iuyi^il^j^i X^(x,).\j{x^)' dXc ' dXc ) 

^ „ (\i(xc)+\i{x^f) ^ / dXiix^) Y I y^ y^ rfA,(xe) _ dXjjXc) 

^i Xi(xc)'^ \ dxc J ^i ^j¥=i dxc dxc 



E l f d\i(xc) \ I v^ f d\i{xc) \ I V^ V^ 

i Xi{x^) ' \ dxc J ^ l^i \ dxc J ^ ^i ^ji- 

_ v^ 1 ( d\i[xc) \ I (s~y 

^i \{x^) \ dxc J ^ \^ 



i dxc dxc 



dXjjXc) 



i dxc 

2 



— \^ 1 f dXi{xc) \ I f _d_ sr^n \ ( ^ \\ _ V^" 1 f dXi(xc) \ f \r)\ 

- l^i A,(X,) ' \ dXr^ ) ^ \dx^ I^i = l^t\.-^C)J - 2^i = l X,(x,) ' \ dXc J y^^) 

where we have used the fact that, for a Poisson distribution, Aj = K{Var{Ii)) = 



-42- 

REFERENCES 

Adorf, H.-M. 1996, Astronomical Data Analysis Software and Systems V, 101, 13 

Anderson, J., & King, I. R. 2000, PASP, 112, 1360 

Auer, L. H., van Altena, W. 1978, AJ, 83, 531 

Bastian, U. 2004. GAIA technical note, 2004BASNOCODE. Available at 

http:www.rssd. esa.intdoc_fetch.php?id=2939027 (last accessed on November 
2012) 

Bendinelli, O., Zavatti, F., & Parmeggiani, G. 1988, Journal of Astrophysics and Astronomy, 
9, 17 

Bristow, P., Piatek, S., & Pryor, C 2005, Space Telescope European Coordinating Facility 
Newsletter, 38, 12 

Bristow, P., Kerber, F., & Rosa, M. R. 2006, The 2005 HST Calibration Workshop: Hubble 
After the Transition to Two-Gyro Mode, 299 

Cramer, H. 1946, Mathematical Methods of Statistics. Princeton University Press, New 
Jersey 

Fruchter, A. S., & Hook, R. N. 2002, PASP, 114, 144 

Gilliland, R. L. 1992, Astronomical CCD Observing and Reduction Techniques, 23, 68 

Howell, S.B. 2006, Handbook of CCD Astronomy, 2nd Edition, Cambridge Observing 
Handbooks for Research Astronomers. Cambridge University Press 

Howell, S.B. 2013. In Astrometry for Astrophysics, William F. van Altena (ed.). Cambridge 
University Press, Chapter 14 



-43- 

Jakobsen, P., Greenfield, P., & Jedrzejewski, R. 1992, A&A, 253, 329 

King, I. R. 1971, PASP, 83, 199 

King, I. R. 1983, PASP, 95, 163 

Lee, J.-F., van Altena, W. 1983, AJ, 88, 1683 

Lindegren, L. 1978, lAU Colloq. 48: Modern Astrometry, 197 

Lindegren, L. 1997, GAIA Technical Note "Astrometric precision for direct fringe detection" 

Lindegren, L. 2005, The Three-Dimensional Universe with Gaia, 576, 29 

Lindegren, L. 2010, ISSI Scientific Reports Series, 9, 279 

Mackay, C. D. 1986, ARA&A, 24, 255 

Mendez, R. A., Costa, E., Pedreros, M. H., et al. 2010, PASP, 122, 853 

Mendez, R. A., Costa, E., Gallart, C, et al. 2011, AJ, 142, 93 

Meyer, S.L. 1992. Data Analysis for Scientists and Engineers, Peer Management Consultants, 
Ltd. 

Monet, D. G. 1992, Astronomical CCD Observing and Reduction Techniques, 23, 221 

Piatek, S., Pryor, C, Olszewski, E. W., et al. 2002, AJ, 124, 3198 

Press, W.H., Flannery B.P., Teukolsky S.A., Vetterling, W.T. 1992. Numerical Recipes in 
Fortran: The Art of Scientific Computing, Cambridge University Press, 2nd edition 

Rao, C. R. 1945, Bulletin of the Calcutta Mathematical Society 37, 81 

Roe, B. P. 2010. Probability and Statistics in Experimental Physics, Springer- Verlag, New 
York, 2nd edition 



-44- 

Stone, R. C. 1989, AJ, 97, 1227 

Stuart, A., Ord, J. K., Arnold S. 2004. Advanced Theory of Statistics, Volume 2A, Classical 
inference and the linear model, Oxford University Press, New York 

van Altena, W. F., & Auer, L. H. 1975, Image Processing Techniques in Astronomy, 54, 411 

Astrometry for Astrophysics, William F. van Altena (ed.). Cambridge University Press 

van Altena, W.F. and Fomalont, E.B. 2013. In Astrometry for Astrophysics, William F. 
van Altena (ed.). Cambridge University Press, Chapter 9 

Zaccheo, T. S., Gonsalves, R. A., Ebstein, S. M., & Nisenson, P. 1995, ApJ, 439, L43 

Winick, K. A., 1986, J. Opt. Soc. Am. A., 3, 1809 



This manuscript was prepared with the A AS lATfrjX macros v5.2. 



-45- 




C^ 



LO 



LO 

d 



09 



Ot 

[SDLU 



02 



- o 







dO 



I 1 

u 

a; 
if) 

o 

D 

I I 

X 



D 



Fig. 1. — Cramer- Rao bound as given be equation ( 12T1) . in milli-arcsec (mas), as a function 
of detector pixel size Ax in arcsec. All curves were computed for a constant background 
per pixel B = 300 ADU/pix and F/B = 10. The dotted, solid, and dashed lines are for a 
Gaussian PSF with a FWHM of 0.5, 1.0, and 1. 5 arcsec r espec tively, all centered in a given 



Winickl ( 119861 )). The dashed-triple dotted 



pixel of width 1.0 pixel (compare to Figure 1 in 

and dashed-dotted lines are for a FWHM of 0.5 arcsec, but de-centered by 0.125 pix and 

0.25 nix from the center of the nixel resnectivelv. 
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Fig. 2. — Cramer- Rao bound as given be equation ( 12T1) . in milli-arcsec (mas), as a function 
of pixel size Ax in arcsec. All curves were computed for a background given by equation fl2^ 
with /, = 2 000 ADU/arcsec, RON = 5e-,D = 0e-,G = 2 e^/ADU, and for a Gaussian 
source with FWHM = 1 arcsec centered on a pixel. The curves shown have different values 
of the flux, and hence a different S/N. From top to bottom we have F = 1 000 ADU, 
S/N ^ 20, F = 2 000 ADU, S/N ^ 35 (both dashed lines); F = 5000 ADU, S/N ^ 74 
fsolid line, as a reference, in this case we have Fmcv ?a Q.SO ADU at Ax = 0.2 arcsec. see 



-47- 



Table 1: Cramer- Rao bound on as a function of pixel offset for a source with FWHM = 
0.5 arcsec observed with two detectors of pixel size indicated by the Ax value. Other pa- 
rameters are those of Fig. [H 



Ax 


Pixel offset 


0"CR 


arcsec 


pix 


mas 


0.5 


0.0 


4.57 


0.5 


0.125 


4.46 


0.5 


0.250 


4.22 


0.7 


0.0 


7.31 


0.7 


0.125 


6.06 


0.7 


0.250 


4.65 


0.9 


0.0 


15.70 


0.9 


0.125 


9.00 


0.9 


0.250 


5.33 
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Table 2: Small-pixel approximation and exact Cramer- Rao bound on some limiting cases. 
Ax FWHM F Fmax /. S/N"" Exact^ Approx^ 

arcsec arcsec ADU ADU ADU mas mas 



0.2 


1.0 


1000 


~186 


2 000 


-20 


27 


24 


0.2 


1.0 


50 000 


-9 308 


2 000 


-300 


1.5 


1.4 


0.2 


1.0 


50 000 


-9 308 


500 


-305 


1.4 


1.4 


0.1 


0.5 


50 000 


-9 308 


4 000 


-300 


0.76 


0.67 


0.1 


0.5 


1000 


- 186 


4 000 


-20 


13 


12 


0.1 


0.5 


5 000 


- 1024 


4 000 


-78 


3.5 


2.4/2. 1<^ 


0.2 


1.0 


5 000 


- 1024 


2 000 


-78 


6.9 


4.8/4.3<^ 



















°As computed from equation [28l with RON = 5 e^, G = 2 e^/ADU, and the tabular value for Aa; and fs- 
'' "Exact" Cramer- Rao bound computed from equation [^ 

■^Small-pixel approximation Cramer- Rao bound computed from equation I45i depending on if S/N is large or 
small. 
■^'The first value is for the case of a weak source, the second for a strong source, both from equation H51 
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Table 3: Comparison of the standard deviation on position a^^ from the ML and LS methods 
and the Cramer- Rao bound, ctcr, as a function of (total) Flux. 









Method: 


ML 


LS 


WLS 


Cramer-Rao 


F 
ADU 


-^max 
ADU 


SjN 




pix 


pix 


pix 


O'CR 

pix 



30, 080 


5600 


-230 


10,002 


1862 


~ 120 


3,222 


600 


-55 


1,612 


300 


-32 


540 


100 


-12 


268 


50 


~ 6 



0.0101 ±0.0005 0.0116 ±0.0006 0.0103 ± 0.0005 0.010 

0.0215 ±0.0010 0.0227 ±0.0010 0.0215 ± 0.0010 0.020 

0.0429 ±0.0019 0.0437 ±0.0020 0.0431 ± 0.0019 0.045 

0.080 ±0.003 0.081 ±0.004 0.080 ± 0.003 0.079 

0.212 ±0.010 0.212 ±0.010 0.212 ± 0.010 0.209 

0.423 ±0.020 0.423 ±0.020 0.423 ± 0.020 0.406 



"The Ict uncertainties in the derived sta ndard devia tion were computed from the variance of the variance, 
Var{u\^), as given by equation (3.44) in lRod (J201Cll ). It is easy to show that the uncertainty in the standard 



deviation would be then given by ^Var[a1.^l2 Gx 



Table 4. Comparison of the best 2-d digital centering algorithms (from IStond (119891 )) and the Cramer- Rao bound as 

a function of (total) Flux. 



FWHM=1 arcsec, B=l ADU/arcsec FWHM=4: arcsec, S=l ADU/arcsec FWHM=1 arcsec, 5=1000 ADU/arcsec FWHM=4. arcsec, 5=1000 ADU/arcsec 
Ax = 0.125 arcsec Ax = 0.5 arcsec Ax = 0.125 arcsec Ax = 0.5 arcsec 



F 


o-st 


0"CR 


o-st 


o"CR 


o-st 


0"CR 


o-st 


o-CR 


ADU 


arcsec 


arcsec 


arcsec 


arcsec 


arcsec 


arcsec 


arcsec 


arcsec 


100 


0.052 


0.047 


0.4 


0.2 





0.24 





1.9 


1,000 


0.015 


0.014 


0.07 


0.055 


0.04 


0.029 


— 


0.20 


10,000 


0.042^ 


0.043 


0.02 


0.017 


0.006 


0.0053 


0.07^^ 


0.027 


100, 000 


0.0012'' 


0.0014 


0.005^ 


0.0054 


0.0015*^ 


0.0014 


0.009== 


0.0060 ^ 



^Moment solution. 

''Far from moment solution, as expected since B is large. 

'^No moment solution found. 



Note. 



All Cramer-Rao estimates used G = 1 c /ADU and RON = e 



